clear all
clc
syms thetai gamma thetaf tauk taun tauc alpha fgdp igdp phi nrate b delta rho


Time=16*330;
taui=1;
%the values that are obtained by averaging the values of the countries in the sample 
alpha=0.4133;
b=0.9658;
delta=0.0461;
    
tauk=0.2457;
taun=0.2263;
tauc=0.1104;
    
fgdp=1638.533;
igdp=352.6533;
phi=1.0671;
nrate=0.2268;
rho=0.025;

labori1=(((1-rho*taui)*gamma*thetai)/((1-taun)*(1-alpha)*thetaf))^(1/(1-gamma));
labori2=(((1/b-1)+delta)/(alpha*(1-tauk)*thetaf))^(alpha/((1-alpha)*(1-gamma)));
labori=labori1*labori2;
laborf1=(Time-labori)*gamma*(1-rho*taui)*thetai*labori^(gamma-1)-phi*(1-rho*taui)*thetai*labori^gamma;
laborf2a=gamma*(1*-rho*taui)*thetai*labori^(gamma-1);
laborf2b=phi*[(alpha*(1-tauk)+(1-alpha)*(1-taun))*thetaf*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(alpha/(1-alpha))-delta*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(1/(1-alpha))];
laborf=laborf1/(laborf2a+laborf2b);
K1=(((1-tauk)*thetaf*alpha)/(1/b-1+delta))^(1/(1-alpha));
K=laborf*K1;
Yf=thetaf*(K^(alpha))*(laborf^(1-alpha));
Yi=thetai*(labori^(gamma));
C=((1-tauk)*alpha*Yf+(1-taun)*(1-alpha)*Yf+(1-rho*taui)*Yi-delta*K)/(1+tauc);
infrate=Yi/Yf;
lrate=labori/laborf;
leisure=Time-laborf-labori;
R=tauc*C+tauk*alpha*Yf+taun*(1-alpha)*Yf+rho*taui*Yi;

%calibration of parameters
z0=[0.5,100,0.5];
options= optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt',...
    'MaxFunctionEvaluations',1500);
res=lsqnonlin(@(z) funcsavg (z,tauk,taun,b,alpha,fgdp,igdp, phi ,nrate,delta,rho),z0, [0,0,0],[],options);
thetaf=res(1);
thetai=res(2);
gamma=res(3);

%% 


clear

Time=16*330;
taui=1;

syms tauk taun tauc thetai gamma thetaf alpha b phi delta rho
labori1=(((1-rho*taui)*gamma*thetai)/((1-taun)*(1-alpha)*thetaf))^(1/(1-gamma));
labori2=(((1/b-1)+delta)/(alpha*(1-tauk)*thetaf))^(alpha/((1-alpha)*(1-gamma)));
labori=labori1*labori2;
laborf1=(Time-labori)*gamma*(1-rho*taui)*thetai*labori^(gamma-1)-phi*(1-rho*taui)*thetai*labori^gamma;
laborf2a=gamma*(1*-rho*taui)*thetai*labori^(gamma-1);
laborf2b=phi*[(alpha*(1-tauk)+(1-alpha)*(1-taun))*thetaf*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(alpha/(1-alpha))-delta*((alpha*(1-tauk)*thetaf)/(1/b-1+delta))^(1/(1-alpha))];
laborf=laborf1/(laborf2a+laborf2b);
K1=(((1-tauk)*thetaf*alpha)/(1/b-1+delta))^(1/(1-alpha));
K=laborf*K1;
Yf=thetaf*(K^(alpha))*(laborf^(1-alpha));
Yi=thetai*(labori^(gamma));



C=((1-tauk)*alpha*Yf+(1-taun)*(1-alpha)*Yf+(1-rho*taui)*Yi-delta*K)/(1+tauc);
R=tauc*C+tauk*alpha*Yf+taun*(1-alpha)*Yf+rho*taui*Yi;
leisure=Time-laborf-labori;
utility=log(C)+phi*log(leisure);
%turevRtaui=diff(R,taui);
%turevUtaui=diff(utility,taui);
turevRtauc=diff(R,tauc);
turevUtauc=diff(utility,tauc);
%mcpftaui=-turevUtaui/turevRtaui;
turevRtauk=diff(R,tauk);
turevUtauk=diff(utility,tauk);
mcpftauk=-(turevUtauk/turevRtauk);
turevRtaun=diff(R,taun);
turevUtaun=diff(utility,taun);
mcpftaun=-(turevUtaun/turevRtaun);
mcpftauc=-(turevUtauc/turevRtauc);



  
    %mcfk=subs(mcpftauk,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    %%mcfc=subs(mcpftauc,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    %mcfn=subs(mcpftaun,{alpha,b,delta,tauk,taun,tauc,phi,thetaf,thetai,gamma},{mcfdata(i,1),mcfdata(i,2),mcfdata(i,3),mcfdata(i,4),mcfdata(i,5),mcfdata(i,6),mcfdata(i,9),mcfdata(i,11),mcfdata(i,12),mcfdata(i,13)});
    alpha=0.4133;
    b=0.9658;
    delta=0.0461;
    
    tauk=0.2457;
    taun=0.2263;
    tauc=0.1104;
    
    fgdp=1638.533;
    igdp=352.6533;
    phi=1.0671;
    nrate=0.2268;
    rho=0.1;
    %values found by calibration
    thetai=24.9994;
   thetaf=0.4325;
    gamma=0.4254;
   mcfdata=zeros(100,11);
   for i=1:100
    rho=i*0.005;   
    mcfdata(i,1)=eval(mcpftauk);
    mcfdata(i,2)=eval(mcpftauc);
    mcfdata(i,3)=eval(mcpftaun);
   
    mcfdata(i,4)=eval(labori);
    mcfdata(i,5)=eval(laborf);
    mcfdata(i,6)=eval(leisure);
    mcfdata(i,7)=eval(K);
    mcfdata(i,8)=eval(C);
    mcfdata(i,9)=eval(R);
    mcfdata(i,10)=eval(Yf);
    mcfdata(i,11)=eval(Yi);
      mcfdata(i,12)=eval(turevUtauk);
    mcfdata(i,13)=eval(turevRtauk);
    mcfdata(i,14)=eval(turevUtaun);
    mcfdata(i,15)=eval(turevRtaun);
    mcfdata(i,16)=eval(turevUtauc);
    mcfdata(i,17)=eval(turevRtauc);
    
    
   end
  % filename = 'avg_outputrho025.xlsx';
%xlswrite(filename,mcfdata);
   
   

